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Abstract 

An epidemic model with distributed time delay is derived to describe the dynamics of in- 
fectious diseases with varying immunity. It is shown that solutions are always positive, and the 
model has at most two steady states: disease-free and endemic. It is proved that the disease- 
free equilibrium is locally and globally asymptotically stable. When an endemic equilibrium 
exists, it is possible to analytically prove its local and global stability using Lyapunov func- 
tionals. Bifurcation analysis is performed using DDE-BIFTOOL and traceDDE to investigate 
different dynamical regimes in the model using numerical continuation for different values of 
system parameters and different integral kernels. 

1 Introduction 

Recent years have witnessed a rapid increase in the use of mathematical models for better under- 
standing of epidemics and disease dynamics. Mathematical models take into account main factors 
that govern development of a disease, such as transmission and recovery rates, and predict how 
the disease will spread over a period of time. Traditional epidemiological models divide the whole 
population into three classes of susceptible, infected and recovered individuals, and the spread of an 
epidemic is governed by the principle of mass action [TJ. For certain diseases, such as, for example, 
sexually transmitted infections, it is important to account for the individuals who have been ex- 
posed to the infection but have not yet become infected, thus the whole population is split into four 
groups, including a separate group of exposed. The incidence rate with which individuals become 
infected is normally taken to be bilinear with respect to the susceptible and infected populations. 
However, there is some evidence that a bilinear incidence rate might not be an effective assumption 
for highly contagious diseases, where a high percentage of the whole population is infected, and the 
transition from susceptible to infected has to be represented by a non-linear function [21 El [3] - 

It is well known that the spread of many infectious diseases can be prevented by vaccination 
of the susceptible population. Furthermore, some infections provide recovered individuals with a 
short or long immunity against re-infection. This means that it is natural to include the effects of 
immunity into the mathematical models in order to better represent the actual dynamics of epidemic 
spread and predict future outbreaks. Immunity can be attained through targeted immunization, 
it can be naturally acquired after an individual has successfully recovered from an infection, and 
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in some cases maternal antibodies can be transmitted to a newborn providing a certain level of 
immunity. In each case the immunity period will vary, as some diseases provide almost life-long 
immunity while others give only a very short-lived non-susceptibility. Quite often, the vaccine- 
induced immunity requires boosting after some period of time, as the vaccine effectiveness wanes 
due to absence of exposure to the disease. In the case of measles, for example, vaccinated individuals 
are less immune than those with naturally acquired immunity [21 E] • Hepatitis B vaccination gives 
only 10 to 15 year immunity and after that the boost is required in order for immunity to remain 
effective [7]. In the case of serogroup C meningococcal disease, the immunity wanes with time 
and its efficacy strongly depends on the immunization programmes available [8] . Recent outbreaks 
of mumps epidemic even amongst vaccinated population suggest that anti-mumps virus antibodies 
wane either due to a wild-type virus or because the population does not have sufficient exposures to 
the disease [9], [10] . Other diseases with waning immunity include varicella virus (although the cause 
of it is still being debated in the literature, some studies suggest that there is a need for a secondary 
vaccination dose |11[ I12j). pertussis, for which immunity declines 6-12 years after the last episode 
of illness or booster dose [T3], and influenza which has a very short lived and strain-dependent 
immunity |14j . 

Delay differential equations have been successfully used to model varying infectious period 
in a range of SIR (susceptible-infected- recovered), SIS (susceptible- infected-susceptible) and SIRS 
(susceptible-infected-recovered-susceptible) epidemic models. Hethcote and van den Driessche have 
considered an SIS epidemic model with variable population size and constant time delay, which 
accounts for duration of infectiousness [15]. They found that an endemic equilibrium, when it 
exists, may undergo Hopf bifurcation and give rise to sustained periodic oscillations for certain 
parameter values. Beretta et al. [16J have studied global stability in an SIR epidemic model with 
distributed delay that describes the time it takes for an individual to lose infectiousness. They have 
used Lyapunov functionals to prove global stability of disease-free and endemic equilibria, provided 
the effective time delay is sufficiently small. Gao et al. [17] have investigated the effects of pulse 
vaccination in an SIR model with distributed delay. They have shown that for sufficiently high 
vaccination rate, the disease-free periodic solution is globally attractive, which implies it is possible 
to completely eradicate the disease from the population. Another interesting integro-differential 
model was studied by Arino et al. [18] . where the authors considered a vaccine, whose effectiveness 
wanes with time according to a general function. They have shown that the system can exhibit 
backward bifurcation and multi-stability for a range of parameter values, which has important 
implication for a design of optimal vaccination policies. A general discussion of time delay models 
in the context of epidemics can be found, for example, in Arino and van den Driessche [19] (and 
references therein). 

Some time ago we put forward a mathematical model of a disease with temporary immunity 
based on a system of delay differential equations [20] . Under quite general assumptions on trans- 
mission incidence, we found conditions for local and global stability of a disease- free and an endemic 
equilibrium when it exists. Several authors have since considered this model with particular choices 
of linear/non-linear incidence rate |21| I22j. and most recently, Brauer et al. have used it to study 
the spread of infection on a two-patch environment [23]. The main assumption was that there is 
a fixed duration of temporary immunity, after which recovered individuals return to the class of 
susceptibles. In this paper, we relax this assumption by considering a more realistic situation, in 
which immunity wanes with time. To model this, we introduce a delay differential equation model 
with distributed delay which takes into account varying immunity period. This means that after 
recovery an individual becomes susceptible again only after spending some time having acquired 
immunity from a disease. As it is important to predict the dynamics of an epidemic, we will prove 
local and global stability of the disease-free equilibrium. Biologically, this stability implies that 
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although the disease may initially be present in the population, as time progresses it will eventually 
die out. Moreover, for some values of parameters it is possible to show that when an endemic 
equilibrium is feasible, it may be globally asymptotically stable. To get a better understanding of 
model dynamics, we will perform numerical bifurcation analysis with a Dirac ^-function kernel us- 
ing a DDE-BIFTOOL package [21]. DDE-BIFTOOL is a numerical continuation tool that detects 
bifurcations and allows one to follow steady states and periodic solutions in the parameter space 
to get a general picture of system's behaviour depending on parameters and, in particular, on the 
time delay. For non-trivial kernels, i.e. weak and strong kernels, we will use traceDDE |25j to find 
the boundary of Hopf bifurcation in terms of system parameters and an average time delay. 

The paper is organised as follows. In Section 2 we introduce a model and prove that solutions of 
the model with positive initial conditions will remain positive for all time. Section 3 is devoted to 
the local and global stability analyses of disease-free and endemic equilibria. Numerical bifurcation 
analysis of the model is performed in Section 4 for several different choices of delay kernels. The 
paper concludes with summary and discussions. 



2 Derivation of the model and positivity of solutions 

In this section we will derive a delayed epidemic model with distributed delay. Recently, Kyrychko 
and Blyuss have introduced and studied a mathematical model for diseases with a constant immu- 
nity time and a nonlinear incidence rate in the form [20j : 



dS ^=V~ <t>mt))S{t) - nS(t) + 7 I(t - T)e-»T, 

d ^ = M(I(t))S{t)-{» + 1 )I{t), (1) 

^ = 7 /(t) - 7 J(t - r)e-^ - »R(t), 

where the total population is divided into three sub-groups, S(t), I(t) and R(t) to denote sus- 
ceptible, infected and recovered individuals, respectively. The parameters are as follows, /i is a 
natural mortality rate, <j) is a disease transmission rate, and 7 is a recovery rate. Nonlinear inci- 
dence rate is represented by the function /(/), and r is a temporary immunity period, after which 
recovered individuals return into the class of susceptible. In this model, the immunity period is 
assumed to be constant and the same for all individuals. It is noteworthy that in the above model 
S(t) + I(t) + R(t) = N(t), and N(t) -> 1 as t -)• oo. 

As it was described in the Introduction, the duration of immunity for different diseases may vary 
from a very short-lived to life-long immunity, and waning rates depend on the amount of exposure, 
boosting times, age etc. To account for this in our model, we assume that immunity period r for 
different individuals varies from to oo. We denote by g(£) the probability density of taking time £ 
to lose acquired immunity, so that g(£)d£ is the probability of losing immunity somewhere between 
£ and £ + dt; after acquiring it. If one normalizes g(£) as follows, 

g(s)ds = 1, and g > 0, (2) 

then the probability of having lost acquired immunity s time units after acquiring it is / g{£)d£, 

Jo 

rs poo 

and the probability of still having immunity s time units after acquiring it is 1 — / g(C)d^ = / <7(£)^£- 

Jo Js 
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The rate of recruitment into the recovered class R at a time t — s is jl(t — s). The probability of 
these individuals still being alive at time t is e~^ s , and the probability of them still being immune 

g(£,)dt;. Integrating over all previous times leads one to the expression 



poo poo pt roc 

R{t) = 7 / I(t- s)e^ s / g(0d£d8 = 7 / I(s)e-^ / g(S)d£d8. (3) 

JO Js J— 00 Jt—s 

Differentiating this expression gives the following integro-differential equation for number of recov- 
ered individuals R(t): 



dR 

Itt 



poo 

jl(t) - 7 / I(t- s)g(s)e-^ s ds - fxR(t). 
Jo 



With this derivation, the model to be investigated in this paper takes the form of a system of delay 
differential equations with distributed delay: 

^ = M - cf>f(I(t))S(t) - nS(t) + 7 I(t - s)g(s)e^ s ds, 

dl{t) 
dt 

dR{t) 
dt 

where, as before, \i is a natural mortality rate, ^> is a disease transmission rate, and 7 is a recovery 
rate. 

It is easy to see that when g(s) = 5(s — r), where 5 is the Dirac delta function, the integral in 
Q becomes time-delayed term in ([I]): 



4>f(I(t))S(t)-(n + y)I(t), (4) 

poo 

-- 7 /(*) - 7 / J(t - s)g(s)e-^ds - »R(t) , 



poo 

/ I{t- s)5{s-T)e-^ds = e-^I{t-T). 
Jo 



In order to simplify system Q we assume that transmission of the disease can be adequately de- 
scribed as bilinear, and, consequently, the term f(I(t))S(t) becomes I(t)S(t). With this assumption 
we arrive at the system 

^ = fi — </>I(t)S(t) - fiS(t) + 7 J™ I{t - s)g(s)e-^ds, 

dl(t) 
dt 

dR(t) 
dt 

where the last equation may be omitted due to the fact that the first two equations do not depend 
on R(t). 

Since the model ^ describes temporal dynamics of human population, it is important to show 
that the solutions of this system do not become negative. We need to show that S(t) > 0, I(t) > 0, 
R(t) > 0, i.e. the solutions of system ([5]) are positive for all t 6 (0;oo). It is more convenient to 



(j>I(t)S(t)-(fx + j)I(t), (5) 

poo 

-- >yl(t) - 7 / I(t- s)g(s)e-^ds - fiR(t), 
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start by proving that I(t) > for all t > 0. Let Iq(s) > 0, s G (— oo; 0), and io(0) > be the initial 
data for I(t). We shall prove the positivity of I(t) by contradiction. Let t\ > be the first time 
when I(t)S(t) = 0. Assuming that I(t{) = implies that S(t) > 0, for all t E [0; t{\. Let 

then, for f E [0;ti], dl/dt > AI(t). Therefore, I{t x ) > I(0)e Atl > 0. This is a contradiction, and, 
hence, I(t) > for all t > 0. Next, let So denote the initial data for S(t), so that S(t) = So(t) for 
all t G (—00; 0). Let So be continuous and satisfy So(i) > for all t G (—00; 0), and <Sb(0) > 0. By 
contradiction, we prove that S(t) > for all t > 0. Assume that there exists a first time to > 
such that S(to) = 0. This indicates that S(t) > for t G [0; to) and 



dt 



/•oo 

/i + 7 y I(t -s) g{s)e~^ s ds > 0. 



>0,Vt e[0;oo) 



Therefore, — ^ °^ > 0. This is a contradiction to our assumption, since it implies S*(t) must be 

negative for t just before 4q, which contradicts the choice of to- The positivity of R(t) for t > 
follows from the integral representation ^ and positivity of I(t). Therefore, we have proved that 
the solutions of system ^ with positive initial conditions remain positive for all times. 



3 Equilibria and their stability 

3.1 Local stability of the disease-free and endemic steady states 



Omitting the last equation in system (|5j), equilibria (S, I) are determined by setting S(t) = I(t) = 0. 
There are two steady states, namely, a disease-free steady state Eq = (1;0) and an endemic steady 
state E = (5,1), where 

S==^and7=- M7 + M)"^ 



7 - [x + 7 / g(s)e ^ s ds 



Using the properties of function g{s) from (J2j), one concludes that while the equilibrium E exists 
for arbitrary values of parameters, it is biologically relevant if and only if 

7 + n < <fi. 

It is important to analyse the stability of these equilibria, as it will indicate whether the disease 
will die out eventually, or it will persist for all time. The linearisation of the system ^ without 
the last equation near the steady state Eq = (1,0) gives 



dS(t) 
dt 

dT(t) 



/•CO 

0J-/xS(t)+7 / T(t- s)g(s)e- lls ds, 

(6) 



^I(t)-(/x + 7)/(t). 



dt 

The characteristic equation for Q has the form 

(A + ^)(A-<A + / u + 7 ) 
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with the corresponding eigenvalues 

Ai = -//; A 2 = 4* ~ M - 7- ( 7 ) 
We define the basic reproduction number as 

v + i 

The basic reproduction number identifies the number of secondary infections from the infected 
source and plays a crucial role in understanding the development of epidemics. From ([7]) it follows 
that when TZq < 1, the disease- free equilibrium Eq = (1,0) is the only steady state of the model, 
and it is locally asymptotically stable. 

When TZo > 1, there exists a non-trivial equilibrium E = (S, I). In order to analyse the stability 
of this equilibrium we use a Lyapunov functional technique. Linearising the system ([5]) without the 
last equation about E by setting S(t) = S + p(t) and I(t) = I + q(t) gives 

dp(t) 



dt 

dqjt) 
dt 



oo 







{-(pi - fi)p(t) - <(>Sq(t) + 7 / q(t - s)g(s)e-' ls ds, 
cj>Ip(t), 



(8) 



where S = — -y — . Introduce the following functional: 

9 

V = -q\t) + -(p(t) + g(t)) 2 + ury / 5 ( S )e- 2 ^ / q 2 (v)dvds, 

where w is a real positive constant, and V(t) > 0. Differentiating V along the solution of system 
Q and substituting q and p we arrive at 

~ />oo 

V= cj)Ipq — wfip 2 — w(n + j)pq + wjp / q(t — s)g(s)e~^ s ds 

Jo 

/•oo 

— wfxpq — w(fi + j)q 2 + wqj / q{t — s)g{s)e~^ s ds 

Jo 

/•oo poo 

+ wiq 2 g(s)e- 2f * s ds-wj g(s) e - 2fis q 2 (t - s)ds. (9) 
Jo Jo 

Estimating the fourth and seventh terms using Cauchy-Schwarz and then Holder inequalities, one 
can rewrite ^ as 



V < pq (pi — wn — w(fi + 7) — wfip 2 — w(j + fi)q 2 H q 2 H p 2 + wjq 2 / g{s)e~ 2 ^ s ds . 

J 2 2 j 



<l 



Since I > 0, we can choose a positive constant w as w = <pl/{2^i + 7). Therefore, the above 
inequality simplifies to 



V < 



7 

""2 



w(p + 9 )• 



As the expression in the right-hand side is negative-definite, provided fi > 7/2, we have proved the 
following theorem: 

Theorem 1. Whenever <f> > fj, + j, the endemic equilibrium E = (S,I) is feasible, and provided 
fj, > 7/2, it is locally asymptotically stable. 
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3.2 Global stability of the disease-free steady state 

It has already been shown that when (j) < fi + 7, the trivial equilibrium Eq of the system ^ is 
locally stable. In fact, in this case we have the following result. 

Theorem 2. // cf> < /i + 7, the disease-free equilibrium Eq = (1,0) is globally asymptotically stable. 

Proof. To prove global stability of the disease-free steady state, we first define N(t) = S(t) + I(t) + 
R(t). Adding all three equations in §5§ gives 

dN 

1F =,-,N(t), 

which implies lim^oo N(t) = 1. Under the assumption <j) < /i + 7, we may choose e > sufficiently 
small, such that (j>{l + e) < fi + 7. Since N(t) — > 1 as t — > 00, for sufficiently large t we have 
S(t) + I(t) + R(t) = N(t) < 1 + e. Then the equation for I(t) becomes 

f t < <t>I(t){l + e-I(t)-R(t)}-(n + 1 )I(t) 
< 0/(t)[l + € -/(t)]-(^ + 7 )/(i) 
= /(t)[0(l + e )-0/(t)-(^ + 7 )]. 

From a simple comparison argument and using the fact that 0(1 + e) < jj, + 7, it follows that 
/(i) — )• as t — t- 00. From equation Q for -R(t) it follows that R(t) — > as t — > 00. Since N(t) —¥ 1, 
it follows that S(t) — > 1 as t — > 00, which implies global asymptotic stability of the trivial steady 
state. ■ 

3.3 Global stability of endemic steady state 

As we have already established, when <fi > /U+7 the disease-free steady state is unstable, and there is 
a feasible endemic equilibrium E = (S, I, R), which according to Theorem 1 is locally asymptotically 
stable. To prove global stability of this steady state, we will employ Lyapunov functional approach. 
Let us introduce new variables ui(t) = S(t) — S, U2(t) = I(t) — I, us(t) = R(t) — R. With this 
change of variables, system ([5]) can be written in the form 

du\ 



js -fjLUi - 4>Su2 - (fiuil + 7 / n 2 (t - s)g(s)e ^ s ds, 



■X) 







-~ = 4>Su2 + (jmxl - (fi + 7)n 2 , (10) 



^=7^2-7/ </•>(/ ->•)//(>•)<• ! ' n<ls - I" 1 -',- 







Global stability of the trivial equilibrium of system (10) implies global stability of an endemic 



equilibrium E of the original system (|5j) . Introduce the following functional 

V(u) = \w(u 1 + U2) 2 + ] 2 {ul + ul), (11) 

where w is a positive constant to be determined later in the calculations. Differentiating V{u) and 
using (fTol) gives 



V(u) = — fiwuf — (w + + 7)^2 — t JbU \ + <j>Su\ 
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+ uiU2[— (n + j)w — fiw + <f)I] + wuij / 1*2 (t — s)g(s)e 113 ds 

Jo 

poo poo 

+ wjU2 / U2(t — s)g(s)e~ >is ds + ^112113 — / 1*2 (t — s)g(s)e~ >is ds. 



It can be easily shown that for non-negative initial conditions S(0) = Sq > 0, I(s) = Io(s) > for 
all s € (-00, 0) with 2(0) = 2 > Ojmd 22(0) = R > 0, one has S(t) < max{l, 50 + 10 + #0} = M 
for alH > [20]. Choosing io = 0// (2[a + 7) in the above equation leads to 

poo 

V < - \iwu\ - [(w + l)Qu + 7) - 4>M]u2 - fml + wu^ / u 2 (t - s)g(s)e~ fis ds 

Jo 

poo poo 

+ W7M2 / i*2 — s)g(s)e~ fMS ds + 7^2^/3 — 7U3 / u>2(t — s)g(s)e~ >MS ds. (12) 



The fourth term can be estimated with the help of Cauchy-Schwartz and Holder inequalities as 
follows 



u>ui7 



U2(t-s)g( S )e^ s d S <^u 2 1 + ^- 



\{t-s)g{s)e- 2 ^ds 



and similar estimates can be made for other muj, i 7^ j terms in (12). This finally gives 



V(u) < 



win 



7\ 2 



M+J ) I) - o.l/ 



«2 



(/i - 7)u 3 + 7 U + 



^(t-s) 5 (s)e- 2 ^ds 



(13) 



Let now the Lyapunov functional for the system ( 10 ) be 



ma) - n») +7 ( w + - 



9{s)e 



-2pLS 



■a\(v)dvds. 



t-s 



Differentiating U and using (13), we obtain 



1 



TJ(u) = V(u) + 1 [w + - )ul{t) / 5 (s)e- 2 ^-7 



w + 



j(t - s)<7(s)e- 2 ^ds 



< -w (// - I) «? - [(/i + I) (w + 1) - 0Af] ul - (/i - 7)^ 



+ 7 w + 



(t) / 5 (s)e- 2 ^ S . 



<i 



Finally, the Lyapunov functional can be estimated as follows 



U (u) < —w ( \i 



7\ „,2 



'//(«; + l)-^-0M 



^2 - (/*-7)«!> 



which is negative-definite if [i > 7 and //(to + 1) — wy/2 — <pM > 0. From Lyapunov-LaSalle type 
theorem (Theorem 2.5.3 of Kuang |26j), it then follows that lim^oo Uk(t)) = 0, k = 1,2,3. Thus, 
we have proved the following result. 

Theorem 3. Let the initial conditions for system |5|) be S(0) = Sq > 0, J(s) = Io(s) > for all 
s e (-oo,0) with 1(0) = 2 > and 22(0) = Rq > 0. 2/ 



/i > 7, (j) > /I + 7, //(w + 1) - ^7/2 - > 0, 



8 




Figure 1: (a) Boundary of Hopf bifurcation of the endemic steady state in terms of r, 7 and fi. 
(b) Two-parameter continuation of Hopf boundary for different values of natural death rate \i. In 
both pictures <p = 10. 

where w = <j>I / (2/z + 7) and M = max{l, 5o + /o + Rq}, then the endemic steady state E of system 
is globally asymptotically stable. 

Since biologically reasonable initial conditions for the system ^ satisfy So + Iq + Rq = 1, the 
conditions of Theorem 3 depend only on system parameters, and the Theorem 3 then holds for all 
initial conditions. 



4 Bifurcations and continuation 

In this section we perform numerical bifurcation analysis of an epidemic model ([5]). The main 
purpose is to get an insight into how the dynamics of the system changes depending on the system 
parameters, and in particular, on the length of the immunity period. We shall proceed by finding 
bifurcation points of the endemic steady state in terms of time delay r, and then continue the 
bifurcating solutions by changing either the transmission rate 4>, the death rate /i, or the recovery 
rate 7. Continuation analysis is an important part in the study of epidemic models, where param- 
eter values carry a lot of uncertainty due the fact that some of them are difficult to obtain from 
experimental data, and some of them vary significantly in the population. Compared to straight- 
forward numerical simulations, bifurcation analysis provides a much more complete picture of the 
underlying dynamical behaviour for a whole range of parameter values. A powerful continuation 
tool for delay differential equations is a MATLAB based tool called DDE-BIFTOOL [M] . It has 
been successfully used for detecting and analysing bifurcations in models of laser dynamics [27] , 
car- following models [28J, neural networks of coupled cells |29j and many other applications. Cur- 
rently, DDE-BIFTOOL does not have capabilities for studying systems with distributed delay, so 
we will use traceDDE package [25] instead. 
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(a) (b) 
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Figure 2: (a) Amplitude of periodic solutions as a function of time delay r. Solid lines denote 
stable periodic orbits, dashed lines denote unstable periodic orbits, (b) Period T of the periodic 
orbits depending on time delay. Parameter values are /i = 0.15 and 7 = 5. 

4.1 Fixed immunity period 

First of all, we consider the case when the integral kernel is the Dirac 5-function g(s) = 5(s — r), 
which corresponds to a fixed period of temporary immunity r. Two particular choices of the 
transmission function /(I) have been analysed, a linear transmission rate f(I) = 7 and a saturated 
transmission rate in the form /(7) = 7/(1 + 7). The results of numerical bifurcation analysis are 
qualitatively similar for both of these functions, hence we only present here the computations for 
/(7) = 7 to complement the analytical results from earlier sections. In all computations, the time 
delay r plays a role of the main bifurcation parameter. After finding an endemic steady state, we 
then detect its bifurcations depending on the values of r and other system parameters. 

Figure [l] shows the boundary of Hopf bifurcation of the endemic steady state E in the space 
of parameters. To find this Hopf boundary, we fix the natural death rate ji and perform a two- 
parameter continuation in terms of a temporary immunity time delay r and a recovery rate 7. It 
is noteworthy that the smaller the death rate fj,, the larger is the region covered by Hopf boundary 
in j-t plane. Furthermore, for each fixed 7, the periodic orbit corresponding to the smaller value 
of r is stable, while its counterpart for the larger value of r is unstable. Below the Hopf boundary, 
the endemic steady state E is stable. 

Next, we take as a starting point a periodic orbit of very small amplitude close to the Hopf 
boundary, and continue this orbit in time delay r. The results of this continuation are shown in 
Fig. [2](a). They indicate that for sufficiently small r, the corresponding periodic orbit is stable 
(solid lines), but then it goes through a fold, and for sufficiently high values of the time delay r 
there is narrow region of a co-existence of stable and unstable periodic orbits (dashed lines). The 
corresponding plot of the periods of those periodic orbits depending on the time delay is shown in 
Fig. [2|b). This figure indicates that the period increases with the time delay, and decreases with 
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Figure 3: Profiles of periodic solutions for \i = 0.15, 7 = 5 and eft = 15. The duration of temporary 
immunity period is (a) r = 1.1739, (b) r = 2.3431, and (c) r = 2.9367. 



the increase of transmission rate <j>. One can explain this observation by the fact that for a higher 
transmission rate (f>, it takes a shorter period of time for individuals to go through a disease cycle 
and return to the class of susceptibles, hence reducing the duration of a periodic cycle. 

The temporary profiles of periodic solutions for different time delays are shown in Fig. [3j where 
all three cases represent stable periodic orbits. For r just above the boundary of Hopf bifurcation, 
the solution shows little variation in either of the variables, but it becomes more pronounced for 
higher r. As is natural to expect, the peak of infectives slightly lags behind the peak of susceptibles, 
and one can also note a certain asymmetry in the profile of the solutions. 



4.2 Gamma distribution 

A more realistic representation of temporary immunity is achieved when the integral kernel is given 
by some non-trivial function g(s). One possible choice is the so-called Gamma distribution delay 
kernel [30] 

9{S)= (n-1)! ' " = 1 > 2 >-> 
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Figure 4: (a) Boundary of the Hopf bifurcation of the endemic steady state with a weak kernel 
9w(s) = ae~ as in terms of r, 7 and [i. The average immunity time is r = 1/a. (b) Two-parameter 
continuation of Hopf boundary for different values of [i. In both pictures (f> = 10. 

where a > is a constant. In this case, one can introduce an average delay as 

poo 

t = I sg(s)ds = n/a. 
Jo 

For simulations, we concentrate on two particular cases: n = 1 and n = 2, which correspond to a 
weak delay kernel and a strong delay kernel, respectively. In the case of weak kernel g w (s) = ae~ as , 
the maximum influx of recovered individuals into the class of susceptibles comes from individuals 
who are currently recovering, while past recoveries have exponentially decreasing contributions. 
On the other hand, for the strong kernel g s (s) = a 2 se~ as , the maximum influx of recovered into 
susceptibles at any time t comes from those who recovered at t — f, where f is the average immunity 
period. 

For numerical bifurcation analysis of system ^ with weak and strong kernels, we use a Matlab 
package traceDDE, which is based on pseudo-spectral differentiation and allows one to find char- 
acteristic roots and stability charts for linear autonomous systems of delay differential equations 
|25j . In all simulations we replaced the upper limit of the integral by a large positive constant L, 
such that e~^ L < 1. 

Figure [4] depicts the boundary of Hopf bifurcation for weak kernel in the space of average time 
delay f , recovery rate 7 and natural mortality [i. Endemic state is stable below the boundary and 
unstable above it. While the dependence on \x is weak, there is a noticeable variation in the critical 
f as a function of 7: the higher is the recovery rate 7, the higher is the average immunity period r, 
at which the endemic steady state loses stability. We note that unlike the case of a Dirac (^-function 
kernel, the dependence of f on 7 at the Hopf boundary is monotonic. 

In Fig. [5] we show Hopf boundary in the case of a strong kernel g s (s)- This figure bears some 
similarity to Figure [l] due to the fact that, as we have already mentioned, the largest contribution 
in the strong kernel comes from the individuals who recovered time f ago, which is reminiscent of 
a (5-function kernel g(s) = 5(s — r). The endemic steady state is stable outside the boundary and 
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Figure 5: (a) Boundary of the Hopf bifurcation of the endemic steady state with a strong kernel 
g s (s) = a 2 se~ as in terms of f, 7 and [i. The average immunity time is f = 2/a. (b) Two-parameter 
continuation of Hopf boundary for different values of [i. In both pictures (ft = 10. 



unstable inside. It is noteworthy that as the natural mortality \x increases, the instability region 
shrinks, and instability occurs for higher values of recovery rate 7. 



5 Conclusions 

In this paper we have derived a delay differential equations model for the dynamics of infectious 
diseases with varying temporary immunity. This model allows temporary immunity to wane with 
time, and is, therefore, a more realistic analogue of previous models which assumed fixed duration of 
immunity period. If transmission coefficient is not too high, the model admits a single disease-free 
equilibrium, which is globally asymptotically stable. For higher values of a disease transmission 
coefficient, there is a feasible endemic equilibrium, whose local and global stability has been proven 
for certain parameter values using a Lyapunov functional approach. 

In order to obtain a better insight into the dynamics of the system in different parameter 
regimes, we have performed numerical bifurcation analysis for several particular choices of the 
integral kernel, including fixed temporary immunity, weak and strong kernels. The results of 
these simulations suggest that endemic equilibrium may lose its stability via a supercritical Hopf 
bifurcation, thus giving rise to stable periodic solutions. Biologically this means that when the 
temporary immunity period is within a certain range, there will be periodic outbreaks of epidemic, 
and the disease will not be eradicated from the population. In the case of delay kernel being a 
Dirac 5-function, as the duration of temporary immunity r increases, these periodic solutions lose 
stability at a fold bifurcation. The period of periodic orbits increases monotonically with r and 
decreases with the increasing disease transmission coefficient. For a weak kernel, the average time 
delay corresponding to the Hopf boundary increases monotonically with the recovery rate 7. The 
case of strong kernel is somewhat similar to the case of Dirac 5-function in that the endemic steady 
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state is unstable for a range of average immunity periods, and the parameter region of instability 
shrinks with the increase of mortality rate /x. 

The model developed in this paper relates to some earlier work on modelling diseases with 
temporary immunity. For example, Cooke and van den Driessche have considered a model with 
constant temporary immunity and latency [31]. They have shown that similarly to our model, 
it is possible to have periodic solutions for some parameter ranges. It is noteworthy that these 
periodic solution occur due to the time delay representing temporary immunity, but there is no 
such effect achieved by considering time delay in latency alone. In comparison, Gomes et al. have 
analysed several ODE models with temporary and partial immunity, as well as vaccination, where 
immunity wanes at a constant rate |32| . They have shown that when a basic reproduction number 
exceeds unity, the solutions either decay linearly to an endemic steady state, or approach it in 
an oscillatory manner. In their findings, immunity waning rate plays an important role in the 
time scale for oscillatory behaviour. While our model also highlights the importance of temporary 
immunity, in contrast to the above model it is also able to sustain periodic solutions describing 
regular epidemic outbreaks. The main feature is that temporary immunity leads to a possible 
destabilization of endemic steady state, and an interesting open question is what effects would 
vaccination have on the dynamics of an epidemic in such situation. 
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